
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <omp.h>
//////////////////////////////
#include <cuda.h>
#include <cuda_runtime.h>
//////////////////////////////

#define REPEAT 256

// Function that runs on GPU
__device__ __inline__ float f(float x,float y){
  return  sinf(cosf(sinf(sqrtf(powf(cosf(x*x+y*y),1.5f)))));
}


// Kernel: apply the function to all elements of an array
__global__ void apply_function(float * __restrict__ a,
			       const float * __restrict__ b){
  int id   = blockIdx.x * blockDim.x + threadIdx.x;
  a[id] = f(a[id],b[id]);
}


int main(void){
  long unsigned int N = 1024L*1024*256;
  float *a = (float *)malloc(N*sizeof(float));
  float *b = (float *)malloc(N*sizeof(float));
  struct timespec ts0,ts1,tsi;
  double time;
  
  fprintf(stderr,"*** GPU ***\n\nFilling arrays\n");
  for (long unsigned int i=0;i<N;i++){
    a[i] = (float)rand()/(float)(RAND_MAX);
    b[i] = (float)rand()/(float)(RAND_MAX);
  }

  clock_gettime(CLOCK_REALTIME,&ts0);
  fprintf(stderr,"Processing arrays\n");

  
  // Allocate arrays on GPU
  float *d_a,*d_b;
  cudaMalloc(&d_a,N*sizeof(float));
  cudaMalloc(&d_b,N*sizeof(float));
  clock_gettime(CLOCK_REALTIME,&tsi);
  time = (tsi.tv_sec-ts0.tv_sec + (tsi.tv_nsec-ts0.tv_nsec)/1e9);
  fprintf(stderr,"\tAllocation       : %.2f\n",time);
  
  // Transfer initial data to GPU
  cudaMemcpy(d_a,a,N*sizeof(float),cudaMemcpyHostToDevice);
  cudaMemcpy(d_b,b,N*sizeof(float),cudaMemcpyHostToDevice);
  clock_gettime(CLOCK_REALTIME,&ts1);
  time = (ts1.tv_sec-tsi.tv_sec + (ts1.tv_nsec-tsi.tv_nsec)/1e9);
  fprintf(stderr,"\tCopy HOST -> GPU : %.2f\n",time);

  // Perform processing on GPU
  for (int count=0;count<REPEAT;count++){
    apply_function<<<N/128,128>>>(d_a,d_b);
  }
  cudaDeviceSynchronize();
  clock_gettime(CLOCK_REALTIME,&tsi);
  time = (tsi.tv_sec-ts1.tv_sec + (tsi.tv_nsec-ts1.tv_nsec)/1e9);
  fprintf(stderr,"\tProcessing       : %.2f\n",time);
 
  // Retrieve the result from GPU
  cudaMemcpy(a,d_a,N*sizeof(float),cudaMemcpyDeviceToHost);
  clock_gettime(CLOCK_REALTIME,&ts1);
  time = (ts1.tv_sec-tsi.tv_sec + (ts1.tv_nsec-tsi.tv_nsec)/1e9);
  fprintf(stderr,"\tCopy GPU -> HOST : %.2f\n",time);

  time = (ts1.tv_sec-ts0.tv_sec + (ts1.tv_nsec-ts0.tv_nsec)/1e9);
  fprintf(stderr,"Total : %.2f seconds\n\n\n",time);

  exit(0);
}
